# === load relevant libraries
library(data.table)
library(ggplot2)
library(DescTools)
library(plm)
library(lfe)
library(dplyr)
library(sandwich)

rm(list = ls())

# academic style portfolios
data = readRDS('input_data/academic_style_data.RDS')
data = data[(yyyymm >= 200201) & (yyyymm <= 200212), list(yyyymm, category, ret, rating)]
data = as.data.table(melt(data, id.vars = c('yyyymm','category')))
names(data)[3] = 'var'

# take out average by period to focus on cross-sectional dispersion
data = merge(data, data[, list(m = mean(value)), list(var, yyyymm)], by = c('yyyymm','var'))
data[, value := value - m]
data[, m := NULL]

# sort by predicted rating changes in 200112
tmp = readRDS('input_data/predicted_academic_style_rating_change_dec_2001.RDS')
data = merge(data, tmp, by = 'category'); rm(tmp)

# enrich for plotting purposes
data[, date := as.Date(as.character(100*yyyymm+1), '%Y%m%d')]
data = data[rank %in% c(1,2,8,9)]
data[, lab := ifelse(rank == 9, 'Most negatively affected',ifelse(rank == 8, 'Second neg',ifelse(rank == 1, 'Most positively affected', 'Second pos')))]
data_bk = copy(data)
event_date = as.Date('2002-06-01')

# output: Figure A7 (a): ratings
data = copy(data_bk[var == 'rating'])
ggplot(data, aes(x = date, y = value, color = reorder(lab, -rank))) + geom_line(lwd = 1) + theme_classic() + 
  labs(x = NULL, y = 'Demeand rating') + theme(legend.position = c(.75, .2), legend.title = element_blank()) + 
  geom_vline(xintercept = event_date, lty = 3) + coord_cartesian(ylim = c(-.8, .6)) + geom_hline(yintercept = 0, lty = 3, lwd = .2)

# output: Figure A7 (b), cumulative returns
data = copy(data_bk[var == 'ret'])
for (this in c(1,2,8,9)){data[rank == this, value := cumsum(log(1+value))]}
rm(this)
ggplot(data, aes(x = date, y = value, color = reorder(lab, -rank))) + geom_line(lwd = 1) + theme_classic() + 
  labs(x = NULL, y = 'Cumulative log return') + theme(legend.position = c(.25, .9), legend.title = element_blank()) + 
  geom_vline(xintercept = event_date, lty = 3) + coord_cartesian(ylim = c(-.15, .25)) + 
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) + geom_hline(yintercept = 0, lty = 3, lwd = .2)
